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Setting upper limits on the strength of periodic gravitational waves 
using the first science data from the GEO 600 and LIGO detectors 
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Data collected by the GEO 600 and LIGO interferometric gravitational wave detectors during 
their first observational science run were searched for continuous gravitational waves from the pulsar 
J1939+2134 at twice its rotation frequency. Two independent analysis methods were used and are 
demonstrated in this paper: a frequency domain method and a time domain method. Both achieve 
consistent null results, placing new upper limits on the strength of the pulsar's gravitational wave 
emission. A model emission mechanism is used to interpret the limits as a constraint on the pulsar's 
equatorial ellipticity. 
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II THE DETECTORS 



I. INTRODUCTION 

This work presents methods to search for periodic 
gravitational waves generated by known pulsars, us- 
ing data collected by interferonietric gravitational wave 
detectors. To illustrate these methods, upper limits 
are placed on the strength of waves emitted by pulsar 
J1939-I-2134 at its expected 1284 Hz emission frequency 
during SI 0. SI is the first observational science run 
of the LIGO HQ and GEO 0,11 detectors and it took 
place during 17 days between August 23 and Septem- 
ber 9, 2002. The sensitivity of the searches presented 
here surpasses that of previous searches for gravitational 
waves from this source. However, measurements of the 
spin-down rate of the pulsar indicate that a detectable 
signal is very unlikely given the instrument performance 
for this data set: for these early observations the detec- 
tors were not operating at their eventual design sensitiv- 
ities. Substantial improvements in detector noise have 
been achieved since the SI observations, and further im- 
provements are planned. We expect that the methods 
presented here will eventually enable the direct detection 
of periodic gravitational waves. 

In Sectional we describe the configuration and calibra- 
tion of the four LIGO and GEO interferometers and de- 
rive their expected sensitivities to periodic sources having 
known locations, frequencies and spindown rates. In Sec- 
tion we consider proposed neutron star gravitational 
wave emission mechanisms, and introduce notation for 
describing the nearly monochromatic signals emitted by 
isolated neutron stars. Statistical properties of the data, 
analyses methods and results are presented in Section lTVI 
These results are then summarized and compared in Sec- 
tion In Section we also interpret the upper limits 
on the signal amplitude as a constraint on the ellipticity 
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of the pulsar and consider our results in the context of 
previous upper limits. 



II. THE DETECTORS 

Gravitational waves are a fundamental consequence 
of Einstein's General Theory of Relativity 0, 0, in 
which they represent perturbations of the spacetime met- 
ric which propagate at the speed of light. Gravita- 
tional waves produced by the acceleration of compact 
astrophysical objects may be detected by monitoring the 
motions they induce on freely-falling test bodies. The 
strength of these waves, called the strain, can be charac- 
terized by the fractional variation in the geodesic sepa- 
ration between these test bodies. 

During the past decade, several scientific collabora- 
tions have constructed a new type of detector for gravita- 
tional waves. These large-scale interferometric detectors 
include the US Laser Interferometer Gravitational Wave 
Observatory (LIGO), located in Hanford, WA, and Liv- 
ingston, LA, built by a Caltech-MIT collaboration 
the GEO 600 detector near Hannover, Germany, built by 
a British-German collaboration 0,0; the VIRGO detec- 
tor in Pisa, Italy, built by an Italian- French collaboration 
[il; and the Japanese TAMA 300 detector in Tokyo 0. In 
these detectors, the relative positions of suspended test 
masses are sensed interferomctrically. A gravitational 
wave produces a time-varying differential displacement 
AL(t) in an interferometer that is proportional to its 
arm length L. The amplitude of the gravitational wave 
is described by the dimensionless strain h{t) = AL{t)/L. 
For realistic periodic astrophysical sources we typically 
expect strain amplitudes smaller than 10~^'^ . 

The following sections introduce the operating config- 
urations of GEO 600 and LIGO detectors during the SI 
run. The references provide more detailed descriptions 
of these detectors. 



A. Instrument configurations 

The GEO 600 detector comprises a 4-beam Michelson 
delay line system of arm length 600 m. The interferom- 
eter is illuminated by frequency-stabilized light from an 
injection-locked Nd:YAG laser. Before reaching the inter- 
ferometer, the light is passed through two 8 m triangular 
mode-cleaning cavities. During SI approximately 2 W of 
light was incident on the interferometer. A power recy- 
cling mirror of 1% transmission was installed to increase 
the effective laser power available for the measurement. 

LIGO comprises three power-recycled Michelson inter- 
ferometers with resonant Fabry-Perot cavity arms. A 
4 km and a 2 km interferometer are collocated at the Han- 
ford site and are designated HI and H2 respectively, and 
a 4 km interferometer at the Livingston site is designated 
LI. Each interferometer employs a Nd:YAG laser sta- 
bilized using a monolithic reference cavity and a 12 m 
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mode-cleaning cavity. 

In all four instruments the beam splitters, recycling 
mirrors and test masses are hung as pendulums from 
multilayer seismic isolation filters to isolate them from 
local forces. The masses and beam paths are housed in 
high vacuum enclosures to preclude optical scintillation 
and acoustic interference. 

Sinusoidal calibration forces of known amplitude were 
applied to the test bodies throughout the observing run. 
These signals were recovered from the data stream and 
used to periodically update the scale factors linking 
recorded signal amplitude to strain. The principal cal- 
ibration uncertainties arise from the imprecision in the 
electromechanical coupling coefficients of the force ac- 
tuators. These were estimated by comparison with the 
known laser wavelength, by actuating a test mass be- 
tween interference fringes. In the Hanford interferome- 
ters, the calibration was also verified against piezoelectric 
displacement transducers connected to the mirror sup- 
port structures. For the SI observations, the net ampli- 
tude uncertainty was estimated at ±4% for GEO, ±10% 
for each of the LIGO interferometers. 



B. Expected sensitivity 

We define the gravitational wave strength ho of a con- 
tinuous signal from a given source as the maximum peak 
amplitude which could be received by an interferome- 
ter if the orientation of the pulsar and the detector were 
both optimal. Thus, ho depends on the intrinsic emission 
strength and on the source distance, but not on the incli- 
nation of the pulsar's spin axis or on the antenna pattern 
of the detector. 

The calibrated interferometer strain output is a time 
series 

s{t) = h{t) + n{t), (2.1) 

where h{t) is the received signal, n{t) is the detector noise 
and t is the time in the detector's frame. 

The noise n(t) is characterized by its single-sided power 
spectral density 5'„(/). Assuming this noise is Gaussian 
and taking some fixed observation time^ T, we can com- 
pute the amplitude of a putative continuous signal 
which would be detected in, e.g., 90% of experimental 
trials if truly present, but would arise randomly from the 
noise background in only 1% of trials (what we call a 1% 
"false alarm rate" and a 10% "false dismissal rate"). 

If we fix a false alarm rate, it is clear that the lower the 
desired false dismissal rate the higher the signal needs to 
be. The detection statistic used in section HV CI provides 
the lowest false dismissal rate for a given false alarm rate 



Here we presume that we know the position, frequency and spin- 
down parameters of the source and that T is between a few days 
and several months. 



and signal strength and it is thus optimal in the Neyman- 
Pearson sense (see for example [13 )■ The amplitude of 
the average signal that we could detect in Gaussian sta- 
tionary noise with a false alarm rate of 1% and a false 
dismissal rate of 10% using the detection statistic de- 
scribed in ^3 is given by^ 

{h„) = 11.4V5„(/,)/r, (2.2) 

where fs is the frequency of the signaP. The upper curves 
in Fig. n show (ho) for the LIGO and GEO detectors 
during SI. Observation times for respective interferom- 
eters are given in the figure. Because of ground motion, 
equipment failures and alignment drifts, the four inter- 
ferometers were not always fully operational during the 
SI run, thus the observation times vary from detector to 
detector. 

The lower curves in Fig.nrcpresent (ho) corresponding 
to the design sensitivity of the various detectors. An 
observation of T = 1 y was assumed. 

The filled circles in Fig. show the constraints that 
measurements of spin-down rates of known pulsars place 
on the expected gravitational wave signal, under the as- 
sumption that the pulsars are rigid rotators with a mo- 
ment of inertia of lO'^^gcm^ and that all the observed 
spin-down rate is due to the emission of gravitational 
waves. 

As shown in Fig. ^ under the above assumptions no 
detection is expected for any known pulsar at the sen- 
sitivity achieved during the SI run. Furthermore, many 
known pulsars are rotating too slowly to be detected by 
the initial ground-based interferometers. However, the 
number of millisecond pulsars observed in this band con- 
tinues to increase with new radio surveys, and the known 
targets plotted here constitute a highly selected sample. 
Future searches for previously undiscovered rotating neu- 
tron stars using the methods presented here will sample 
a different and potentially much larger subset of the total 
population. 



III. PERIODIC GRAVITATIONAL WAVES 

A. Expected emission by neutron stars 

The strongest argument that some neutron stars (NSs) 
are emitting gravitational waves (GWs) with amplitude 



^ The average is over different positions, inchnations and polariza- 
tions of the source^ 

^ This differs from Il2ll for three reasons: 1) the ho used here is 
twice that defined in |12|| . 2) we use a different statistic for this 
detection problem (a chi-square distribution with 4 degrees of 
freedom) and 3) we have spe cified a false dismissal rate of 10% 
whereas the derivation in |12| has an implicit false dismissal rate 
of about 50%. If we use this false dismissal rate and the 
statistic we get (ho) = 7.61/ Sn{fs)/T. 
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frequency (Hz) 



FIG. 1: Upper curves: characteristic amplitude (ho) of a 
known monochromatic signal detectable with a 1% false alarm 
rate and a 10% false dismissal rate by the GEO and LIGO 
detectors at the SI sensitivity and with an observation time 
equal to the up-time of the detectors during SI (GEO: 401 h, 
Ll:137h, Hl:209h, H2:214h). Lower curves: {ho) for the 
design sensitivities of the detectors for an assumed 1 y obser- 
vation time. Filled circles: upper limit on (ho) from measured 
spin-down rate of known radio pulsars assuming a moment of 
inertia of 10*^ g cm^. These upper limits were derived under 
the assumption that all the measured loss of angular momen- 
tum of the star is due to the emission of gravitational waves, 
neglecting spm-down contribution from electromagnetic and 
particle emission. The arrow points to the filled circle repre- 
senting pulsar J1939-I-2134. 

detectable by Advanced LIGO 0, /lo > lO^^^ _ 19-26^ 
is due to Bildsten 0,^3- He noted that the inferred ro- 
tation frequencies of low-mass X-ray binaries (LMXBs) 
are all clustered in the range fr ~ 270 — 620 Hz (an infer- 
ence strengthed by the recent observations of flil It^). 
whereas a priori there should be no cut-off in /r, up to the 
(estimated) NS break-up frequency of ~ 1.5 kHz. Updat- 
ing a suggestion by Wagoner [l^ [13 , Bildsten proposed 
that LMXBs have reached an equilibrium where spin- 
up due to accretion is balanced by spin-down from GW 
emission. Since the GW spin-down torque scales like f^, 
a wide range of accretion rates then leads to a rather 
narrow range of equilibrium rotation rates, as observed. 

Millisecond pulsars (MSPs) are generally believed to 
be recycled pulsars: old pulsars that were spun up by ac- 
cretion during an LMXB phase fj20l f2l|). The rotation 
rates of MSPs also show a high-frequency cut-off ; the 
fastest (PSR J1939+2134) has = 642 Hz. If the GWs 
that arrest the spin-up of accreting NSs continue to be 
emitted in the MSP phase (e.g., because of some per- 
sistent deformation of the NS shape away from axisym- 
mctry), then they could also account for the observed 
spin-down of MSPs. In this case, the GW amplitudes 



of MSPs would in fact be (very close to) the 'spin-down 
upper limits' shown in Fig.^ (Note that the MSP's spin- 
down rate is generally attributed entirely to the pulsar's 
magnetic field; indeed, pulsar magnetic fields are typi- 
cally inferred this way. However, there appears to be no 
strong evidence supporting this inference.) 

We now turn to the possible physical mechanisms 
responsible for periodic GWs in this frequency range. 
The main possibilities that have been considered are 1) 
NS spin precession, 2) an excited NS oscillation mode 
(mostly likely the r-mode), and 3) small distortions of 
the NS shape away from axisymmetry. At present the 
third mechanism (small ellipticity) seems the most plau- 
sible source of detectable GWs, and in this paper we 
set upper limits for this particular mechanism (the three 
mechanisms predict three different GW frequencies for 
the same observed rotation frequency). However, we be- 
gin by briefly commenting on the other two possibilities. 

A NS precesses (or 'wobbles') when its angular mo- 
mentum J is not aligned with any principal axis of its 
inertia tensor. A wobbling NS emits GWs at the inertial- 
frame precession frequency, which is very nearly the rota- 
tion frequency, f^-. While large-amplitude wobble could 
plausibly produce GW amplitudes ho ^ 10~^^ over short 
timescales, the problem with this mechanism is that dis- 
sipation should damp NS wobble quickly while this 
dissipation timescale is quite uncertain (it is perhaps of 
order a year for a MSP), it is almost certainly orders of 
magnitude shorter than typical lifetimes of MSPs. So 
unless some mechanism is found that regularly re-excites 
large-amplitude wobble, it is unlikely that any nearby 
MSP would be wobbling. Moreover, most MSPs have 
highly stable pulse shapes, and typically appear not to 
be wobbling substantially. In particular, the single-pulse 
characteristics of PSR J1939-I-2134 have been observed 
to be extremely stable with no pulse-to-pulse variation 
except for occasional giant pulses |2^. It has been ver- 
ified through radio observations that PSR J1939-I-2134 
continued to spin according to a simple spin-down model 
during SI [H. 

R-modes (modes driven by Coriolis forces) have been 
a source of excitement among GW theorists since 1998, 
when Andersson and Friedman and Morsink 
showed that they should be unstable due to gravitational 
back-reaction (the Chandrasekhar-Friedman-Schutz in- 
stability). Nonlinear mode- mode coupling is predicted 
to saturate the growth of r-modes at dimensionlcss am- 
plitude a < lO^^if./kRz)^^^ [23. This implies r-mode 
radiation from nascent NSs in extragalactic supernovae 
will not be detectable, but r-mode GWs from old, re- 
cycled Galactic NSs could still be detectable by Ad- 
vanced LIGO. For example, GWs from an excited r-mode 
could balance accretion torque in accreting NSs, as in the 
Wagoner-Bildsten mechanism. 

We now turn to GWs from small non-axisymmetries in 
the NS shape. If hg is the amplitude of the signal at the 
detector from an optimally oriented source, as described 
above, and if we assume that the emission mechanism is 
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due to deviations of the pulsar's shape from perfect axial 
symmetry, then 



n-o — e, 



(3.1) 



where r is the distance to the NS, Izz is its principal 
moment of inertia about the rotation axis, e = {I^x — 
Iyy)/Izz is its ellipticity, and the gravitational wave signal 
frequency, fs, is exactly twice the rotation frequency, /,-. 
G is Newton's constant, and c is the speed of light. This 
is the emission mechanism that we assume produces the 
gravitational wave signal that we are targeting. 

One possible source of ellipticity is tiny 'hills' in the 
NS crust, which are supported by crustal shear stresses. 
In this case, the maximum ellipticity is set by the crustal 
breaking strain (T,„ax : 



5x 10-8(a^,,/10- 



(3.2) 



The coefficient in Ea. l3.2l is low both because a NS's crust 
is rather thin (compared to the NS radius), and because 
the crust's shear modulus fi is small compared to the 
ambient pressure p: p/p ~ 10^'^ — 10^^. (If NSs have 
solid cores, as well as crusts, then much larger elliptic- 
ities would be possible.) For the LMXBs, Ushomirsky, 
Cutler & Bildsten showed that lateral temperature 
variations in the crust of order 5%, or lateral composi- 
tion variations of order 0.5% (in the charge-to-mass ra- 
tio), could build up NS ellipticities of order 10~® — 10"^, 
but only if the crust's breaking strain is large enough to 
sustain such hills. 

Strong internal magnetic fields are another possible 
source of NS ellipticity. Cutler has argued that if a 
NS's interior magnetic field B has a toroidal topology (as 
expected if the B field was generated by strong differen- 
tial rotation immediately after collapse) , then dissipation 
tends to re-orient the symmetry axis of the toroidal B 
field perpendicular to the rotation axis, which is the ideal 
orientation for maximizing equatorial ellipticity. Toroidal 
B fields of order lOi^-lQi^ G would lead to sufficient GW 
emission to halt the spin-up of LMXBs and account for 
the observed spin-down of MSPs. 



B. The signal received from an isolated pulsar 

A gravitational wave signal we detect from an isolated 
pulsar will be amplitude-modulated by the varying sen- 
sitivity of the detector as it rotates with the Earth (the 
detector's 'antenna pattern'). The detected strain has 
the form 

, / N „ , N , 1 + COS^ L , , 

h{t) = F+(t, ^) /lo ^ cos$(t) 

+ F^{t,ij) ho cost sin$(t), (3.3) 

where t is the angle between neutron star's spin direction 
s and the propagation direction of the waves k , and (f>(t) 



right ascension (J2000) ig*" 39"" 38".560 210(2) 
declination (J2000) -^21° 34'° 59M41 66(6) 
RA proper motion —0.130(8) masyr"'^ 
dec proper motion —0.464(9) masyr^^ 

period (1//0 0.001 557 806 468 819 794(2) s 
period derivative 1.051 193(2) x 10"^^ ss"^ 
epocli of period and position MJDN 47 500 



TABLE I: Parameters for the target pulsar of the analy- 
ses presented here, PSR J1939-f2134 (also designated PSR 
B1937-f21). Numbers in parentheses indicate uncertainty in 
the last digit. 



is the phase evolution of the signal, i^+.x sltc the strain 
antenna patterns of the detector to the plus and cross 
polarizations and are bounded between -1 and 1. They 
depend on the orientation of the detector and the source, 
and on the polarization of the waves, described by the 
polarization angle ip ^. 

The signal will also be Doppler shifted by the orbital 
motion and rotation of the Earth. The resulting phase 
evolution of the received signal can be described by a 
truncated Taylor series as 



$(t) = ^o + 2tt[MT-To) 



^fs{T-To) 



IfsiT- 



(3.4) 



where 



T = t + St = t 



rd • k 



Abo - A 



so- 



(3.5) 



Here, T is the time of arrival of a signal at the solar sys- 
tem barycentre (SSB), (po is the phase of the signal at 
fiducial time Tq, rd is the position of the detector with 
respect to the SSB, and A_eo and As© are the solar sys- 
tem Einstein and Shapiro time delays respectively [30|. 

The timing routines used to compute the conversion 
between terrestrial and SSB time (Eq. I3.5|l were checked 
by comparison with the widely-used radio astronomy 
timing package TEMPO 0|. This comparison (Figure|3 
confirmed an accuracy of better than ± 4 ps, thus ensur- 
ing no more than 0.01 radian phase mismatch between a 
putative signal and its template. This results in a negli- 
gible fractional signal-to-noise ratio loss, of order ~ 10~^. 

Table |2 shows the parameters of the pulsar that we 
have chosen to illustrate our analysis methods [s^ . 



* Following the conventions of llll , ip is the angle (clockwise about 
k) from z X k to k X s, where z is directed to the North Celestial 
Pole, k X s is the i-axis of the wave frame - also called the waves' 
principal + polarization direction. 
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10= 



8x10" 



4x10'' 



2x10* 




spin-down 



-8.663 3(43) x lO^^Hzs" 



-10 



-5 5 

timing residuals (yiis) 



FIG. 2: Histogram of timing residuals between our bary cen- 
tring routines and TEMPO, derived by comparing the phase 
evolution of test signals produced by the two software pack- 
ages. 156 locations in the sky were chosen at random and 
the residuals calculated once an hour for the entire year 2002. 
The maximum timing error is < 4 /is. 



IV. DATA ANALYSES 

A. Introduction 

Two independent search methods are presented here: 
i) a frequency domain method which can be employed for 
exploring large parameter space volumes and ii) a time 
domain method for targeted searches of systems with an 
arbitrary but known phase evolution. 

Both approaches will be used to cast an upper limit on 
the amplitude of the periodic gravitational wave signal: 
a Bayesian approach for the time domain analysis and a 
frequentist approach for the frequency domain analysis. 
These approaches provide answers to two different ques- 
tions and therefore should not be expected to result in 
the exact same numerical answer |33 . ^3 • The frequen- 
tist upper limit refers to the reliability of a procedure for 
identifying an interval that contains the true value of hg. 
In particular, the frequentist confidence level is the frac- 
tion of putative observations in which, in the presence of 
a signal at the level of the upper limit value identified 
by the actual measurement, /ig^^", the upper limit identi- 
fied by the frequentist procedure would have been higher 
than /iq^^°. The Bayesian upper limit, on the other hand, 
defines an interval in ho that, based on the observation 
made and on prior beliefs, includes the true value with 
95% probability. The probability that we associate with 
the Bayesian upper limit characterizes the uncertainty in 
ho given the observation made. This is distinct from the 
reliability, evaluated over an ensemble of observations, of 
a procedure for identifying intervals. 

All the software used for the analyses is part of the 
publicly available LSC Algorithm Library (LAL, [13 )• 
This is a library that comprises roughly 700 functions 
specific to gravitational wave data analysis. 



/s at start of GEO observation 1 283.856 487 705(5) Hz 

/s at start of LI observation 1 283.856 487 692(5) Hz 

/s at start of HI observation 1 283.856 487687(5) Hz 

/s at start of H2 observation 1 283.856 487682(5) Hz 



TABLE H: Run parameters for PSR J1939-I-2134. The dif- 
ferent emission frequencies correspond to the different initial 
epochs at which each of the searches began. Numbers in 
parentheses indicate the uncertainty in the last digit or digits. 



B. Statistical characterization of the data 

Due to the narrow frequency band in which the tar- 
get signal has appreciable energy, it is most convenient 
to characterize the noise in the frequency domain. We 
divided the data into 60 s blocks and took the Fourier 
transform of each. The resulting set of Fourier trans- 
forms will be referred to as SFTs (Short time-baseline 
Fourier Transforms) and is described in more detail in 

HvTTTI 

The frequency of the pulsar signal at the beginning of 
the observation for every detector is reported in Table ITU 
Also reported is the value of the spin-down parameter 
expressed in units of Hzs~^. We have studied the statis- 
tical properties of the data in a narrow frequency band 
(0.5 Hz) containing the emission frequency. This is the 
frequency search region, as well as the region used for 
estimating both the noise background and the detection 
efficiency. Fig. |31 summarizes our findings. Two types 
of distributions are plotted. The first column shows the 
distributions of bin power; for each SFT (labelled by a) 
and for every frequency bin (labelled by 1 < fc < M) in 
the band 1283.75 to 1284.25 Hz, we have computed the 
quantity 



Pak = 



\Xak\ 



(4.1) 



where Xak is the SFT datum at frequency index k of 
the a-th SFT, and have histogrammed these values. If 
the data are Gaussian and if the different frequency bins 
in every SFT are independent realizations of the same 
random process, then we expect the normalized power 
variable described above (Pak) to follow an exponential 
distribution with a mean and standard deviation of 1, 
as shown by the dashed line. The circles are the experi- 
mental points. The standard deviation of the measured 
distribution for GEO data is 0.95. The LIGO Livingston, 
Hanford 4 km and Hanford 2 km data are also shown in 
Fig. 13 The standard deviation of the Pak for all of these 
is 0.97. 

The plots in the second column of Fig.|21show the dis- 
tribution of phase differences between adjacent frequency 
bins. With the same notation as above, we have com- 
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FIG. 3: Histograms of Pak and of A^ak for the four detectors. 



puted the quantity 



fe-i, 



(4.2) 



where is the phase of the SFT datum at frequency 
index k of the a-th SFT and the difference is reduced to 
the range [— tt, tt]. A$afc is therefore the distance in phase 
between data at adjacent frequency bins. If the data were 
from a purely random process we expect this distribution 
to be uniform between — tt and tt, as observed. 

Fig.01shows the average value of y/S^ over a 1 Hz band 
from 1283.5 to 1284.5 Hz as a function of time in days 
for the entire SI run starting from the beginning of SI 
(15:00 UTC, August 23 2002). These plots monitor the 
stationarity of the noise in the band of interest over the 
course of the run. 

Fig. shows v^SVi as a fimction of frequency between 
1281 and 1 285 Hz. During SI the received signal is ex- 
pected to have a frequency 1 283.8 Hz. This frequency 
is shown as a dashed vertical line. During the SI obser- 
vation time the Doppler modulation changed this signal 
frequency by no more than 0.03 Hz, two SFT frequency 
bins. For these plots 5'„ has been estimated by averaging 
the power in each frequency bin over the entire SI run. 
A broad spectral feature is observed in the GEO data. 
This feature is 0.5 Hz wide, comparatively broad with 
respect to the expected Doppler shift of the target sig- 
nal, and represents only a 10% perturbation in the local 
power spectral density. 




10 15 20 
days 



days 



FIG. 4: The square root of the average value of Sn for all 
four interferometers over a band of 1 Hz starting at 1 283.5 Hz 
versus time in days starting at the beginning of SI (August 
23 2002, 15:00 UTC). 
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FIG. 5: in a band of 4 Hz (starting at 1281Hz) using 

the entire SI data set analyzed from the four interferometers. 
The noise is shown in units of 10"^° Hz"^/^. The dashed 
vertical line indicates the expected frequency of the signal 
received from J1939-I-2134. 



C. The frequency domain technique 

1. The Short time baseline Fourier Transforms 

In principle the only constraint on the time baseline of 
the SFTs used in the frequency domain analysis is that 
the instantaneous frequency of a putative signal does not 
shift during the time baseline by more than half a fre- 
quency bin. For frequencies in the kiloHertz range this 



implies a maximum time baseline of order 30min (hav- 
ing assumed an observation time of several months and a 
source declination roughly the same as the latitude of the 
detector). However, in practice, since we are also estimat- 
ing the noise on the same time baseline, it is advisable 
for the time baseline to be short enough to follow the 
non-stationarities of the system. On the other hand, for 
the frequency-domain analysis, the computational time 
required to carry out a search increases linearly with the 
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number of Fourier transforms. Thus the shorter the time 
basehne, the higher the computational load. We have 
chosen for the SI run a time baseline of 60s as a com- 
promise between the two opposing needs. 

Interruptions in interferometer operation broke each 
time series into segments separated by gaps representing 
invalid or contaminated data. Only valid data segments 
were included in the analysis. Each valid 60s data seg- 
ment was filtered with a fifth-order Butterworth high- 
pass filter having a knee frequency of 100 Hz. Then, a 
nearly flat-top Tukey window function was applied to 
each data segment in the time domain. The window 
changes the value of less than 1% of the data in each 
60s chunk. Each data segment was then Fast Fourier 
Transformed and written to an SET file. These SETs 
were computed once and then used repeatedly for differ- 
ent analyses. 



2. The T statistic 

The detection statistic that we use is described in ■ 
As in we call this statistic JF^, though differences 
between our definition and that given in are pointed 
out below. 

The T statistic derives from the method of maximum 
likelihood. The log likelihood function, In A, is, for Gaus- 
sian noise 

In K^{s\h)-\{h\h), (4.3) 

where 

W,) = «|"£(|md/. ,4.4) 

s is the calibrated detector output time series, h is the 
target signal (commonly referred to as the template) , is 
the Fourier transform operator, and 5'„(/) is the one- 
sided power spectral density of the noise. The T statis- 
tic is the maximum value of In A with respect to all the 
unknown signals parameters, given our data and a set of 
known template parameters. In fact, if some or all of the 
signal's parameters are unknown, it is standard practice 
to compute the likelihood for different template parame- 
ters and look for the highest values. The maximum of the 
likelihood function is the statistic of choice for matched 
filtering methods, and it is the optimal detection statis- 
tic as defined by the Neyman-Pearson criterion: lowest 
false dismissal rate at a fixed false alarm rate (see, for 
example. Section III B|) . 



^ Note that this detection statistic has nothing to do with the F- 
statistic of the statistical literature, which is ratio of two sample 
variances, or the F-test of the null hypothesis that the two sam- 
ples are drawn from distributions of the same variance. 



In our case the known parameters are the position of 
the source (a, 5 angles on the celestial sphere), the emis- 
sion frequency /s and the first order spin-down parameter 
value /g. The unknown parameters are the orientation of 
the pulsar (angle /,), the polarization state of the wave 
(angle ?/;), its initial phase (/)o, and the wave's amplitude 

The core of the calculation of J- consists in computing 
integrals of the type given in Eq. 14.41 using templates for 
the two polarizations of the wave. The results are opti- 
mally combined as described in |lll | except we consider a 
single frequency component signal. Also, we do not treat 
Snif) as constant in time. Thus, while the method is de- 
fined in |ll| in the context of stationary Gaussian noise, 
we adapt it so that it can be used even when the noise is 
nonstationary. The calculation is easily performed in the 
frequency domain since signal energy is concentrated in 
a narrow frequency band. Using the SETs described in 
II V G II some approximations can be made to simplify the 
calculation and improve computational efficiency while 
still recovering most (> 98%) of the signal power. 

The method of computing J- was developed for a spe- 
cific computational architecture: a cost-effective Beowulf 
cluster, which is an ensemble of loosely-coupled proces- 
sors with simple network architecture. This becomes 
crucial when exploring very large parameter-space vol- 
umes for unknown sources using long observation peri- 
ods, because the search depth and breadth are limited 
by computational resources. The SI analyses described 
here were carried out using Gondor [33 | on the Merlin 
and Medusa clusters at the AEI and UWM respectively 
[Ullil- Each cluster has 300 independent GPUs. 

As a point of reference we note that it takes of order 
of a few seconds of GPU time on a 1.8 GHz-class GPU 
to determine the statistic for a single template with 
16 d of observation time. 



3. Setting an upper limit on ho 

The outcome JF* of a specific targeted search repre- 
sents the optimal detection statistic for that search. Over 
an independent ensemble of similar searches in the pres- 
ence of stationary Gaussian noise, 2J-* is a random vari- 
able that follows a distribution with four degrees of 
freedom. If the data also contain a signal, this distribu- 
tion has a non-centrality parameter A proportional to the 
time-integral of the squared signal. 

Detection of that signal would be signified by a large 
value J-* unlikely to have arisen from the noise-only dis- 
tribution. If instead the value is consistent with pure 
noise (as we find in this instance), we can place an upper 
limit on the strength of any signal present, as follows: 

Let J^* be the value of the detection statistic in our 
actual experiment. Had there existed in the data a real 
signal with amplitude greater than or equal to /io(C), 
then in an ensemble of identical experiments with dif- 
ferent realizations of the noise, some fraction of trials C 
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would yield a detection statistic exceeding the value J-*. 
We will therefore say that we have placed an upper limit 
ho (C) on the strength of the targeted signal, with confi- 
dence C. This is a standard frequentist upper limit. 

To determine the probability distribution p{2J-\ho), we 
produce a set of simulated artificial signals with fixed 
amplitude ho from fictional pulsars at the position of our 
target source, and with the same spin-down parameter 
value, but with intrinsic emission frequencies that differ 
from it by a few tenths of a hertz. We inject each of 
these artificial signals into our data and run a search 
with a perfectly matched template. For each artificial 
signal we obtain an independent value of the detection 
statistic; wc then histogram these values. If the SFT 
data in nearby frequency bins (of order 100 bins) can be 
considered as different realizations of the same random 
process (justified in llV B|) , then it is reasonable to assume 
that the normalized histogram represents the probability 
density function p{2J-'\ho). One can then compute the 
confidence 

/•OO 

C{ho)= / p{2T\ho)d{2T), (4.5) 

where ho(C) is the functional inverse of C{ho). In prac- 
tice, the value of the integral in Eq. 14.51 is calculated 
directly from our simulations as follows: we count how 
many values of are greater or equal to J-* and divide 
this number by the total number of values. The value 
derived in this way does not rely on any assumptions 
about the shape of the probability distribution function 
(pdf) curve p(2jr|/io). 

There is one more subtlety that must be addressed: 
all eight signal parameters must be specified for each in- 
jected artificial signal. The values of source position and 
spin-down parameters are known from radio data and 
are used for these injections. Every injected signal has 
a different frequency but all such frequencies lie in bins 
that are close to the expected frequency of the target sig- 
nal, 1 283.86 Hz. The values l and ip are not known, and 
no attempt has been made in this analysis to give them 
informative priors based on radio data. However, the 
value of the non-ccntrality parameter that determines the 
p{2J-\hQ) distribution does depend on these values. This 
means that for a given a different confidence level 
can be assigned for the same signal strength, depending 
on the choice of l and "0. 

There are two ways to proceed: either inject a popula- 
tion of signals with different values of t and ip, distributed 
according to the priors on these parameters^, or pick a 
single value for l and for ip. In the latter case it is rea- 
sonable to choose the most pessimistic orientation and 
polarization of the pulsar with respect to the detector 
during the observation time. For fixed signal strength, 



this choice results in the lowest confidence level and thus, 
at fixed confidence, in the most conservative upper limit 
on the signal strength. We have decided to use in our sig- 
nal injection the worst-case values for l (which is always 
7r/2) and ^p, i.e., the values for which the non-centrality 
parameter is the smallest. 

4- The frequency domain SI analysis for PSR J1939+2134 

Table HhI summarizes the results of the frequency do- 
main analysis. For every interferometer (column 1) 
the value of the detection statistic for the search for 
J1939-I-2134 is reported: 2!F*, shown in column 4. Next 
to it is the corresponding value of the chance probability: 

/•OO 

Po{2T*)= / p{2T\ho = 0)d{2T), (4.6) 

our estimate of how frequently one would expect to ob- 
serve the measured value of J-* or greater in the absence 
of a signal. The values of Po{2J^*) are not significant; we 
therefore conclude that there is no evidence of a signal 
and proceed to set an upper limit. 

Tobs is the length of the live-observation time, /iq"^'''^* is 
the amplitude of the population of injected signals that 
yielded a 95% confidence. The upper limit differs 
from /i™''*''^* only by the calibration uncertainty, as ex- 
plained in Section IIV El Cmcas is the confidence level 
derived from the injections of fake signals, and AC its 
estimated uncertainty due to the finite sample size of the 
simulation. 

The quantities in the remaining columns can be used to 
evaluate how far the reported results are from those that 
one expects. The results shown are remarkably consistent 
with what one expects based on the noise and on the 
injected signal: the confidence levels that we determine 
differ from the expected ones by less than 2%. 

Given a perfectly matched template, the expected non- 
centrality parameter when a signal h{t) is added to white 
noise with spectral density Sn is 



where U = Jj, ^ \h{t)\'^ dt. U can also be computed by 
feeding the analysis pipeline pure signal and by perform- 
ing the search with a perfectly matched template ^ hav- 
ing set Sn{f) = 1 s. In Table Hill we report the values of 
Uo, for the worst-case h{t) signals for PSR J1939-F2134 
as 'seen' by the interferometers during their respective 
observation times and with Hq = 2 x 10~^^. The dif- 
ferent values of Uq refiect the different durations of the 



This is indeed one of the consistency checks that have been per- 
® The time domain analysis assumes uniform priors on cos t and formed to vaUdate the analysis software. We have verified that 

i/). the two values of U agree within a 1% accuracy. 
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IFO 


Tobs [d] 


.inject 




Po(2J^*) 


{l/S„)-i [Hz-i] 


;7o/10-»^ [s] 


-^cxp 


-^best-fit 


C'cxp 


Cbest-fit 


C'mcas dz AC 


GEO 


16.7 


1.94 X 10"^^ 


1.5 


0.83 


5.3 X 10"^* 


1.0 


3.6 


3.3 


95.7% 


95.2% 


95.01 ± 0.23% 


LI 


5.73 


2.70 X 10~^^ 


3.6 


0.46 


1.4 X 10"^" 


0.37 


9.6 


8.3 


96.7% 


95.0% 


95.00 ± 0.23% 


HI 


8.73 


5.37 X 10"^^ 


6.0 


0.20 


5.4 X 10"""' 


0.5 


13.3 


12.8 


96.6% 


95.0% 


95.00 ± 0.23% 


H2 


8.90 


3.97 X 10"^^ 


3.4 


0.49 


3.8 X 10"""' 


0.45 


9.3 


7.9 


96.8% 


95.0% 


95.00 ± 0.23% 



TABLE III: Summary of the frequency domain search analyses. Tobs indicates the total duration of the analyzed data set. J^* 
is the measured value of the detection statistic. Po{2J^*) is the probability of getting this value or greater by chance, i.e. in 
the absence of any signal. /Iq"^'^'^' is the amplitude of the population of fake signals that were injected in the data set such that, 
when searched for with a perfectly matched template, Cmcas% of the time the resulting value of was greater than JT*. {l/5„) 
is the average value of the inverse of the noise in a small frequency band around the target frequency. Uo is the time integral 
of the square of the targeted signal with an amplitude of 2 x lO"'^^, at the output of the interferometers, for observations 
times equal to Tobs and in the absence of noise. Aexp is the value of the non-centrality parameter that one expects for the 
distribution of from searches with perfectly matched templates on a population of injected signals with amplitude /ig"-''"^' 
and noise with average power (1/Sn)~^- Abest-flt is the best-fit non-centrality parameter value derived from the distribution 
p{2J^\h''^^'"^^) derived from the software signal injections and searches with perfectly matched templates. Cexp and Cbest-fit are 
the corresponding confidence values for J^* . 



observations and the different orientations of each detec- 
tor with respect to the source. The expected value of the 
non-centrality parameter can be estimated as: 

(.inject \ ^ 
2TiF^) ■ ^^-^^ 

If the noise were stationary, then S„ may be easily de- 
termined. Our noise is not completely stationary, so the 
value determined for the non-centrality parameter A is 
sensitive to the details of how Sn is estimated. The value 
of (l/Sn) used to determine the expected value of A is 
computed as 

I 



where the frequency index k varies over a band ^ 0.2 Hz 
around 1283.89 Hz. N and At are the number of sam- 
ples and the sampling time of the 60s time series that 
are Fourier transformed. Wc choose an harmonic mean 
rather than an arithmetic mean because this is the way 
Sn enters the actual numerical calculation of the J-- 
statistic. This method is advantageous because the es- 
timate it produces is relatively insensitive to very large 
outliers that would otherwise bias the estimate. 

Aoxp is the expected value of the non-centrality param- 
eter based on Sn and /iq"^""^*. Abest-flt is the best-fit value 
of the non-centrality parameter based on the measured 
distribution of values from the simulation. Cexp and 
Cbest-fit are the confidence levels corresponding to these 
distributions integrated between 2J-* and oo. 



Fig. El shows the distributions for p{2T\hQ^^'"^^) . The 
circles result from the simulations described above. The 
solid lines show the best fit non-central curves. The 
shaded region is the integral ofp{2T\h^Q^'"^^) between 2T* 
and oo. By definition this area is 0.95. 



D. The time domain search technique 

1. Overview 

Frequency domain methods offer high search efficien- 
cies when the frequency of the signal and/or the position 
of the neutron star are unknown and need to be deter- 
mined along with the other signal parameters. However 
in the case of known pulsars, where both the intrinsic ro- 
tation frequency of the neutron star and its position are 
known to high accuracy, alternative time domain meth- 



ods become attractive. At some level the two domains are 
of course equivalent, but issues such as data dropouts and 
the handling of signals with complicated phase evolutions 
can be conceptually (and practically) more straightfor- 
ward in a time series analysis than in an analysis based 
on Fourier transforms. 

The time domain search technique employed here in- 
volves multiplying (heterodyning) the quasi-sinusoidal 
signal from the pulsar with a unit-amplitude complex 
function that has a phase evolution equal but of oppo- 
site sign to that of the signal. By carefully modelling this 
expected phase, $(i), we can take account of both the in- 
trinsic frequency and spin-down rate of the neutron star 
and its Doppler shift. In this way the time-dependence 
of the signal is reduced to that of the strain antenna pat- 
tern, and we are left with a relatively simple model-fitting 
problem to infer the unknown pulsar parameters ho, l, ip 
and 00 defined in Eas 13.31 and 13.41 
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FIG. 6: Measured pdf for 2J^ for all four interferometer data 
with injected signals as described in Table IHll The circles 
represent the measured pdf values from the Monte Carlo sim- 
ulations. The lines represent distributions with 4 degrees 
of freedom and best-fit non centrality parameters given in 
Table UTTI The filled area represents the integral of the pdfs 
between 2J^* and -|-cx3. 



In the time domain analysis we take a Bayesian ap- 
proach and therefore express our results in terms of pos- 
terior probability distribution functions for the parame- 
ters of interest. Such pdfs are conceptually very different 
from those used to describe the T statistic used in the 
frequency domain search, and represent the distribution 
of our degree of belief in the values of the unknown pa- 
rameters, based on the experiments and stated prior pdfs. 

The time domain search algorithm comprises stages 
of heterodyning, noise estimation and parameter esti- 
mation. In outline, the data are first heterodyned at 
a constant frequency close to the expected frequency of 
the signal, low-pass filtered to suppress contamination 
from strong signals elsewhere in the detector band and re- 
binned to reduce the sampling frequency from 16 384 Hz 
to 4 Hz. A second (fine) heterodyne is applied to the data 
to account for the time-varying Doppler shift and spin- 
down of the pulsar and any final instrumental calibration, 
and the data are re-binned to 1 sample per minute. We 
take the data as stationary during this period, and make 
an estimate of the noise variance in each 1 min bin from 
the variance and covariancc of the data contributing to 
that bin. This variance is used in the likelihood function 
described below. 

The parameter estimation stage, at which we set the 
Bayesian upper limit on Hq, proceeds from the joint prob- 
ability of these 1-min complex samples, {Bk}- We take 
these Bk values to have a Gaussian likelihood with re- 
spect to our signal model, y{tk] a), where a is a vector in 
our parameter space with components (Hq, (j)o) and 
tk is the time-stamp of the fcth sample. The signal model. 



the complex heterodyne of Eq. 13.31 is 



-Fx(ife;V)/io cos6e''2'^°. 



(4.10) 



We choose uniform prior probabilities for 0o over [0, 27r] 
and "0 over [— 7r/4, 7r/4], and a prior for t that is uniform 
in cosi over [—1,1], corresponding to a uniform proba- 
bility per unit solid angle of pulsar orientation. These 
uniform priors are uninformative in the sense that they 
arc invariant under changes of origin for the parameters. 
Although strictly a scale parameter, the prior for Hq is 
also chosen as constant for Hq > 0, and zero Hq < 0. This 
is a highly informative prior, in the sense that it states 
that the prior probability that ho lies between 10"^'' and 
10"^^ is 10 times less than the prior probability it lies be- 
tween 10^^"^ and 10^^"*, but guarantees that our posterior 
pdf can be normalized. 

The joint posterior pdf for these parameters is 



p(a|{Sfe}) cx p(a)cxp 



E 



SR{gfc-y(tfc;a)}^ 



X exp 



2at 



(4.11) 



where p(a) (cx sint) is the prior on a, cr^^^^j. is the 
variance of the real parts of Bk, and cr^^^^y is the vari- 
ance of the imaginary parts of Bk- 

The final stage in the analysis is to integrate this pos- 
terior pdf over the l, and (po parameters to give a 
marginalized posterior for ho of 



p{ho\{Bk})(x JJJ p{a\{Bk})didi;d^o, (4.12) 

normalized so that p{ho\{Bk}) dho = 1. This curve 
represents the distribution of our degree of belief in any 
particular value of /iq, given the model of the pulsar sig- 
nal, our priors for the pulsar parameters, and the data. 
The width of the curve roughly indicates the range in 
values consistent with our state of knowledge. 

By definition, given our data and priors, there is a 
probability of 0.95 that the true value of ho lies below 
where 



0.95 = 



p{ho\{Bk})dho, 



(4.13) 



and this defines our 95%-credible Bayesian upper limit 
on ho- 

An attraction of this analysis is that data from dif- 
ferent detectors can be combined directly using the ap- 
propriate signal model for each. The combined posterior 
distribution from all the available interferometers comes 
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naturally out of a Bayesian analysis and, for independent 
observations, is simply the (normalized) product of the 
contributing probability distributions, i.e., 

p(a|all data) oc p(a|GEO) x p(a|Hl) x p(a|H2) x p(a|Ll). 

(4.14) 

This posterior pdf embodies all we believe we know about 
the values of the parameters, optimally combining the 
data from all the interferometers in a coherent way. For 
interferometers with very different sensitivities, this will 
closely approximate the result from the most sensitive 
instrument. Again, we must marginalize over t, ip and (pQ 
to obtain the posterior pdf for Hq alone. We note that this 
is more than simply a combination of the marginalized 
pdfs from the separate interferometers as the coherence 
between the instruments is preserved, and it recognizes 
the different polarization sensitivities of each. 

Equipment timing errors discovered after SI cautioned 
against a coherent multi-interferometer analysis with this 
data set. In principle we could assign a suitable prior for 
the resulting phase offsets and marginalize over them. 
However, the dominant position of the Livingston 4km 
interferometer means that even a fully a coherent analysis 
would only improve our sensitivity by about 20%, so we 
have not pursued this. Fully coherent analyses will be 
possible in future observing runs. 



2. The time- domain SI analyses for PSR J1939+2134 

The time domain search used contiguous data seg- 
ments 300 s or longer in duration. 

The effectiveness of the noise estimation procedure de- 
scribed above was assessed from histograms of B/a = 

^{Bk)/o'3t{Bk} +'^{Bk)/<^'3{Bk}- If ^^"^ estimates are cor- 
rect, and our likelihood function is well modelled by a 
Gaussian, these histograms (Fig. [Tj) should also be Gaus- 
sian with a variance of one. Since we divide the noise 
between the real and imaginary components we expect 
the value of to be close (within V2N) of the num- 
ber of real and imaginary data, N (twice the number 
of complex binned data values Bk). A small number of 
outliers with magnitudes of B^ j larger than 5 were not 
included in this or subsequent analyses. 

The marginalized posterior pdfs for /iq are plotted as 
the solid lines in Fig.|51 These represent the distribution 
of our degree of belief in the value of ft-o: following SI, for 
each interferometer. The width of each curve roughly in- 
dicates the range in values consistent with our priors and 
the data from the instruments individually. The formal 
95% upper limits from this analysis are the upper bounds 
to the shaded regions in the plots, and are 2.2 x 10~^^ 
for GEO, 1.4 x lO'^^ for LI, 3.3 x lO'^^ for HI, and 
2.4 X 10-22 for H2. 

The dotted hue in the GEO plot of Fig. El shows the 
(very different) marginalized posterior pdf obtained when 
a simulated signal is added to the data with an amplitude 
of 2.2 X 10-2\ and with = 0°, V = 0° and i = 0°. 
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FIG. 7: Histograms of B/cr = 3fJ(Bfe)/(Jjf{i3,}-f 9(Bfc)/'^a{Bfc} 
for each interferometer. The dotted lines represent the ex- 
pected Gaussian distribution, with = and (7 = 1. 




FIG. 8: For each interferometer, the solid line represents the 
marginalized posterior pdf for /lo (PSR J1939-I-2134) resulting 
from the SI data. The 95% upper limits (extent of the shaded 
region) are 2.2 x 10"^^ for GEO, 1.4 x 10"^^ for LI, 3.3 x 
10"^^ for HI and 2.4 x 10"^^ for H2. The dotted line in the 
GEO plot shows the posterior pdf of /lo in the presence of a 
simulated signal injected into the GEO SI data stream using 
/lo = 2.2 X 10-^\ 00 = 0°, V = 0° and t = 0°. 



Here there is a clear non-zero lower limit for the value 
of /iQ, and a result such as this would have indicated a 
nominal detection, had we seen it. 
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E. Estimation of uncertainties 

In the frequency domain analysis the uncertainty on 
the upper hmit value, /iq^^, has two contributions. The 
first stems from the uncertainty on the confidence (AC « 
0.23%) that results from the finite sample size of the sim- 
ulations. In order to convert this uncertainty into an un- 
certainty on /ig^'^", we have performed several additional 
Monte Carlo simulations. For every run we have injected 
a population of signals with a given strength, /iq"^""^*, near 
/ig^'^° , searched for each of them with a perfectly matched 
template and derived a value of T. With these values we 
were able to estimate the Hq^C) curve near /ig^^", its slope 
/ig and from this the uncertainty on the value of 

Ahf^" w h'oAC. (4.15) 

The second contribution to the uncertainty on the 
value of /ig^*^" comes from errors in the calibration of 
the instruments, which influence the absolute sensitiv- 
ity scale. In particular this reflects in an uncertainty in 
the actual value of the strength of injected signals so that 
^95% _ /j^J'"^* -|- r^Yic sum of this error, estimated 

m lira and the error arising from the finite sample size, 
Eq. 14.151 is given in the frequentist results in Table IIVI 

Note that when a pulsar signal is present in the data, 
errors in the calibration introduce errors in the phase and 
amplitude of that signal. The errors in due to the signal 
are quadratic with the errors in the phase and are linear 
with the errors in the amplitude. However the estimate 
of the noise spectral density is also effected by calibration 
errors, and in particular by the amplitude errors. The net 
effect on J- is that the resulting error on this quantity 
(which can be considered a sort of signal-to-noise ratio) 
is quadratic in calibration errors, thus insensitive, to first 
order, to calibration errors. 

The errors quoted for the Bayesian results in Table Hvl 
simply reflect the calibration uncertainties given in III Al 
For clarity, no attempt has been made to fold a prior for 
this calibration factor into the marginal analysis. 

V. CONCLUSION 
A. Summary of results 

Table HVl summarizes the 95% upper limit (UL) results 
that we have presented in the previous sections. We 
should stress once more that the two analyses address 
two well-posed but different questions, and the common 
nomenclature is somewhat misleading. 

The frequentist upper limit statements made in Section 
IIV CI refer to the likelihood of measuring a given value of 
the detection statistic or greater in repeated experiments, 
assuming a value for ho and a least-favorable orientation 
for the pulsar. The Bayesian limits set in Section llV D II 
refer to the cumulative probability of the value of ho itself 
given the data and prior beliefs in the parameter values. 



V CONCLUSION 



IFO 


Frequentist FDS 


Bayesian TDS 


GEO 


(1.9 ±0.1) X 10"^^ 


(2.2 ±0.1) X 10" 


21 


LI 


(2.7 ±0.3) X 10~^^ 


(1.4 ±0.1) X 10" 


22 


HI 


(5.4 ±0.6) X 10"^^ 


(3.3 ±0.3) X 10" 


22 


H2 


(4.0 ±0.5) X 10"^^ 


(2.4 ±0.2) X 10" 


22 



TABLE IV: Summary of the 95% upper limit values of ho 
for PSR J1939±2134. The frequency domain search (FDS) 
quotes a conservative frequentist upper limit and the time 
domain search (TDS) a Bayesian upper limit after marginal- 
izing over the unknown t, and 00 parameters. 



The Bayesian upper limits report intervals in which we 
are 95% certain that the true value resides. We do not 
expect two such distinct definitions of 'upper limit' to 
yield the same numerical value. 

Recall that the frequentist UL is conservative: it is 
calculated for the worst-case values of signal parameters 
L and "ip ■ The Bayesian TDS method marginalizes over 
these parameters, gathering together the evidence sup- 
porting a particular /ig irrespective of orientation. We 
have also performed an alternative calculation of the fre- 
quentist ULs by using a p{T\ho) derived from a popu- 
lation of signals with cos t and tp parameters uniformly 
distributed, as were the Bayesian priors in the time do- 
main search. As expected, we find that the resulting ULs 
have somewhat lower values than the conservative ones 
reported in TablcHVl 1.2x10^21^ 1.5 x lO'^^, 4.5 x 10"22 
and 2.3 x 10"22 for the GEO, LI, HI and H2 data sets 
respectively. 



B. Discussion of previous upper limit results 

Two prior upper limits have been published on the 
strain of a signal from our specific pulsar J1939+2134. A 
limit oih< 3.1 x lO-^"^ and 1.5 x lO-^"^ for the first and 
second harmonic of the rotation frequency of the pulsar, 
respectively, was set in using 4d of data from the 
Caltech 40 m interferometer. A tighter limit h < 10"^° 
was determined using a divided-bar gravitational wave 
detector at Glasgow University for the second harmonic 
alone 0- 

More sensitive untargeted UL results on the strain 
of periodic GW signals at other frequencies come from 
acoustic bar detector experiments , and . Due 
to the narrow sensitivity bands of these detectors (less 
than 1 Hz around each mode) , and the fact that their fre- 
quencies do not correspond to those of any known pul- 
sars®, studies with bar antennas have not investigated 



* With the exception of the Austrahan detector NIOBE and of the 
Japanese torsional antenna built specifically to detect periodic 
signals from the Crab pulsar 1431 . 
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possible emission from any known pulsars. 

In a UL of 2.9 x 10"^'* was reported for periodic 
signals from the Galactic center, with 921.32 < fs < 
921.38 Hz and no appreciable spin-down over ^ 95.7 days 
of observation. These data were collected by the EX- 
PLORER detector in 1991. This UL result was not ob- 
tained by a coherent search over the entire observation 
time, due to insufficient timing accuracy. 

In ^3 a fully coherent 2 day-long all-sky search was 
performed again on 1991 EXPLORER data in a fs search 
band of about 1 Hz centered at 922 Hz and including 1 
spin-down parameter. It resulted in an UL of 2.8 x 10~^^ 
at the 99% confidence level. This search was based on 
the same detection statistic used in our frequency domain 
analysis. 

Another parameter space search is presented in |42l| . 
Data taken from the ALLEGRO detector during the first 
three months of 1994 were searched for periodic gravita- 
tional wave signals from the Galactic center and from 
the globular cluster 47Tuc, with no resolvable spin-down 
and with fs in the two sensitive bands of their antenna, 
896.30- 897.30Hz and 919.76 - 920.76Hz, with a 10 ^Hz 
resolution. The resulting UL at 8 x lO"^"* is reported. 

There exist several results from searches using e arly 
broadband interferometric detectors 0, 0, EH EM 
Due to the poor sensitivities of these early detector 
prototypes, none of these upper limits is competitive with 
the strain sensitivity achieved here. However, many of 
the new issues and complications associated with broad- 
band search instruments were first confronted in these 
early papers, laying the foundations for future analyses. 

Data from the first science run of the TAMA detector 
were searched for continuous waves from SN1987A in a 
0.05 Hz band at ~ 934.9 Hz. The reported 99% confi- 
dence upper limit was ~ 5 x 10"^'^ |50l |. 

Improved noise performance and longer observation 
times achieved with interferometric detectors since SI has 
made their sensitivities comparable to or better than the 
narrow band peak sensitivity of the acoustic bars cited 
above, over much broader bandwidths. Combined with 
the advances in analysis methods presented in this paper 
we anticipate significant advances in search depth and 
breadth in the next set of observations. 



spindown rate: e < 3.80 x 10-^(10'^^ gcm^/ 1 ^z Y/^. How- 
ever, an ellipticity of ^ 10~^ could in principle be gener- 
ated by an interior magnetic field of strength ~ 10^^ G, 
or it could probably be sustained in a NS with a solid 
core. Therefore, the above exercise suggests that with 
improved detector sensitivities, even a null result from 
a search for unknown pulsars will place interesting con- 
straints on the ellipticities of rapidly-rotating neutron 
stars that might exist in our galactic neighborhood. 
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C. Upper limit on the ellipticity of the pulsar 

An UL on ho for J1939-f 2134 can be interpreted as an 
UL on the neutron star's equatorial ellipticity. Taking 
the distance to J1939-I-2134 to be 3.6 kpc. Eg. 10 elves 
an UL on cUipticity corresponding to /iq^'^" < 1.4 x 10~^^ 
of 

e^^^" = 2.9 X 10- ^l^^^llJ^^ . (5.1) 

Of course, the UL on the cUipticity of J1939+2134 
derived from SI data is about five orders of magnitude 
higher than the UL obtained from the pulsar's measured 
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